home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / integration / tests.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-05-05  |  5.9 KB  |  221 lines
  1. /* integration/tests.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <math.h>
  22. #include <gsl/gsl_math.h>
  23.  
  24. #include "tests.h"
  25.  
  26. /* These are the test functions from table 4.1 of the QUADPACK book */
  27.  
  28. /* f1(x) = x^alpha * log(1/x) */
  29. /* integ(f1,x,0,1) = 1/(alpha + 1)^2 */
  30.  
  31. double f1 (double x, void * params) {
  32.   double alpha = *(double *) params ;
  33.   return pow(x,alpha) * log(1/x) ;
  34. }
  35.  
  36. /* f2(x) = 4^-alpha / ((x-pi/4)^2 + 16^-alpha) */
  37. /* integ(f2,x,0,1) = arctan((4-pi)4^(alpha-1)) + arctan(pi 4^(alpha-1)) */
  38.  
  39. double f2 (double x, void * params) {
  40.   double alpha = *(double *) params ;
  41.   return pow(4.0,-alpha) / (pow((x-M_PI/4.0),2.0) + pow(16.0,-alpha)) ;
  42. }
  43.  
  44. /* f3(x) = cos(2^alpha * sin(x)) */
  45. /* integ(f3,x,0,pi) = pi J_0(2^alpha) */
  46.  
  47. double f3 (double x, void * params) {
  48.   double alpha = *(double *) params ;
  49.   return cos(pow(2.0,alpha) * sin(x)) ;
  50. }
  51.  
  52. /* Functions 4, 5 and 6 are duplicates of functions  1, 2 and 3 */
  53. /* ....                                                         */
  54.  
  55. /* f7(x) = |x - 1/3|^alpha */
  56. /* integ(f7,x,0,1) = ((2/3)^(alpha+1) + (1/3)^(alpha+1))/(alpha + 1) */
  57.  
  58. double f7 (double x, void * params) {
  59.   double alpha = *(double *) params ;
  60.   return pow(fabs(x - (1.0/3.0)),alpha) ;
  61. }
  62.  
  63. /* f8(x) = |x - pi/4|^alpha */
  64. /* integ(f8,x,0,1) = 
  65.    ((1 - pi/4)^(alpha+1) + (pi/4)^(alpha+1))/(alpha + 1) */
  66.  
  67. double f8 (double x, void * params) {
  68.   double alpha = *(double *) params ;
  69.   return pow(fabs(x - (M_PI/4.0)),alpha) ;
  70. }
  71.  
  72. /* f9(x) = sqrt(1 - x^2) / (x + 1 + 2^-alpha) */
  73. /* integ(f9,x,-1,1) = pi/sqrt((1+2^-alpha)^2-1) */
  74.  
  75. double f9 (double x, void * params) {
  76.   double alpha = *(double *) params ;
  77.   return 1 / ((x + 1 + pow(2.0,-alpha)) * sqrt(1-x*x)) ;
  78. }
  79.  
  80. /* f10(x) = sin(x)^(alpha - 1) */
  81. /* integ(f10,x,0,pi/2) = 2^(alpha-2) ((Gamma(alpha/2))^2)/Gamma(alpha) */
  82.  
  83. double f10 (double x, void * params) {
  84.   double alpha = *(double *) params ;
  85.   return pow(sin(x), alpha-1) ;
  86. }
  87.  
  88. /* f11(x) = log(1/x)^(alpha - 1) */
  89. /* integ(f11,x,0,1) = Gamma(alpha) */
  90.  
  91. double f11 (double x, void * params) {
  92.   double alpha = *(double *) params ;
  93.   return pow(log(1/x), alpha-1) ;
  94. }
  95.  
  96. /* f12(x) = exp(20*(x-1)) * sin(2^alpha * x) */
  97. /* integ(f12,x,0,1) = 
  98.    (20 sin(2^alpha) - 2^alpha cos(2^alpha) + 2^alpha exp(-20))
  99.    /(400 + 4^alpha) */
  100.  
  101. double f12 (double x, void * params) {
  102.   double alpha = *(double *) params ;
  103.   return exp(20*(x-1)) * sin(pow(2.0,alpha) * x) ;
  104. }
  105.  
  106. /* f13(x) = cos(2^alpha * x)/sqrt(x(1 - x)) */
  107. /* integ(f13,x,0,1) = pi cos(2^(alpha-1)) J_0(2^(alpha-1))  */
  108.  
  109. double f13 (double x, void * params) {
  110.   double alpha = *(double *) params ;
  111.   return cos(pow(2.0,alpha)*x)/sqrt(x*(1-x)) ;
  112. }
  113.  
  114. double f14 (double x, void * params) {
  115.   double alpha = *(double *) params ;
  116.   return exp(-pow(2.0,-alpha)*x)*cos(x)/sqrt(x) ;
  117. }
  118.  
  119. double f15 (double x, void * params) {
  120.   double alpha = *(double *) params ;
  121.   return x*x * exp(-pow(2.0,-alpha)*x) ;
  122. }
  123.  
  124. double f16 (double x, void * params) {
  125.   double alpha = *(double *) params ;
  126.   if (x==0 && alpha == 1) return 1 ;  /* make the function continuous in x */
  127.   if (x==0 && alpha > 1) return 0 ;   /* avoid problems with pow(0,1) */
  128.   return pow(x,alpha-1)/pow((1+10*x),2.0) ;
  129. }
  130.  
  131. double f17 (double x, void * params) {
  132.   double alpha = *(double *) params ;
  133.   return pow(2.0,-alpha)/(((x-1)*(x-1)+pow(4.0,-alpha))*(x-2)) ;
  134. }
  135.  
  136. /* f454(x) = x^3 log|(x^2-1)(x^2-2)| */
  137. /* integ(f454,x,0,inf) = 61 log(2) + (77/4) log(7) - 27 */
  138.  
  139. double f454 (double x, void * params) {
  140.   double x2 = x * x;
  141.   double x3 = x * x2;
  142.   params = 0 ;
  143.   return x3 * log(fabs((x2 - 1.0) * (x2 - 2.0))) ;
  144. }
  145.  
  146. /* f455(x) = log(x)/(1+100*x^2) */
  147. /* integ(f455,x,0,inf) = -log(10)/20 */
  148.  
  149. double f455 (double x, void * params) {
  150.   params = 0 ;
  151.   return log(x) / (1.0 + 100.0 * x * x) ;
  152. }
  153.  
  154. /* f456(x) = log(x) */
  155. /* integ(f456*sin(10 pi x),x,0,1) = -(gamma + log(10pi) - Ci(10pi))/(10pi) */
  156.  
  157. double f456 (double x, void * params) {
  158.   params = 0 ;
  159.   if (x == 0.0)
  160.     {
  161.       return 0;
  162.     }
  163.   return log(x) ;
  164. }
  165.  
  166. /* f457(x) = 1/sqrt(x) */
  167. /* integ(f457*cos(pi x / 2),x,0,+inf) = 1 */
  168.  
  169. double f457 (double x, void * params) {
  170.   params = 0 ;
  171.   if (x == 0.0)
  172.     {
  173.       return 0;
  174.     }
  175.   return 1/sqrt(x) ;
  176. }
  177.  
  178. /* f458(x) = 1/(1 + log(x)^2)^2 */
  179. /* integ(log(x) f458(x),x,0,1) = (Ci(1) sin(1) + (pi/2 - Si(1)) cos(1))/pi 
  180.                                = -0.1892752 */
  181.  
  182. double f458 (double x, void * params) {
  183.   params = 0 ;
  184.  
  185.   if (x == 0.0) 
  186.     {
  187.       return 0;
  188.     }
  189.   else 
  190.     {
  191.       double u = log(x);
  192.       double v = 1 + u * u;
  193.       
  194.       return 1.0 / (v * v) ;
  195.     }
  196. }
  197.  
  198. /* f459(x) = 1/(5 x^3 + 6) */
  199. /* integ(f459/(x-0),x,-1,5) = log(125/631)/18 */
  200.  
  201. double f459 (double x, void * params) {
  202.   params = 0 ;
  203.   return 1.0 / (5.0 * x * x * x + 6.0) ;
  204. }
  205.  
  206. /* myfn1(x) = exp(-x - x^2) */
  207. /* integ(myfn1,x,-inf,inf) = sqrt(pi) exp(-1/4) */
  208.  
  209. double myfn1 (double x, void * params) {
  210.   params = 0;
  211.   return exp(-x - x*x) ;
  212. }
  213.  
  214. /* myfn2(x) = exp(alpha*x) */
  215. /* integ(myfn2,x,-inf,b) = exp(alpha*b)/alpha */
  216.  
  217. double myfn2 (double x, void * params) {
  218.   double alpha = *(double *) params ;
  219.   return exp(alpha*x) ;
  220. }
  221.